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I.   INTRODUCTION 

This  thesis  presents  and  discusses  a  computer  program, 
called  N ATI -3D,  which  employs  the  finite  element  method  to 
set  up  and  solve  the  elasto dynamic  equations  for  the  free 
vibration  of  a  three-dimensional,  linearly  elastic  body  so 
as  to  emerge  with  an  accurate  approximation  of  the  funda- 
mental frequency  of  free  harmonic  vibration. 

The  investigation  was  originally  motivated  by  the 
interesting,  important,  and  unsolved  problem  of  the  effect 
of  foundation  compliance  and  foundation  inertia  upon  the 
fundamental  natural  frequency  of  a  short  elastic  cantilever 
beam.   It  is  clear  that  the  usual  assumption  that  the  root 
of  such  a  beam  is  rigidly  fixed  does  not  adequately  represent 
the  fact  that  no  foundation  is  truly  rigid.   Calculations 
based  upon  ideal  root  fixity  thus  necessarily  result  in 
estimates  of  the  fundamental  vibration  frequency  that  are 
too  high.   There  seems  to  be  no  reasonably  reliable  way  of 
estimating  this  eff ect  [Ref . 1 ] . 

During  the  course  of  the  work  involved  in  adapting 
existing  computer  software  to  the  problem  indicated  in  the 
last  paragraph  (in  particular  adding  a  dynamic  capability 
which  such  software  did  not  previously  possess)  it  became 
more  and  more  apparent  that  it  would  be  most  useful  to 
focus  attention  upon  the  development  of  a  general  purpose 
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.  program  for  the  evaluation  of  the  fundamental  frequency  of 
three  dimensional  elastic  bodies  of  any  shape.   Because  of 
the  usual  and  not  unexpected  problems  of  testing  and  debugging 
a  very  complicated  digital  computer  program,  time  was  not 
available  for  the  study  of  the  root  fixity  problem.   Further- 
more, and  more  importantly,  for  this  problem  or  any  problem 
in  which  truly  complicated  geometry  is  involved,  digital 
computer  programs  which  rely  solely  on  in-core  storage,  as 
does  the  program  developed  herein  as  well  as  most  of  its 
predecessors,  are  capable  of  supplying  results  of  marginal 
accuracy  so  that  parametric  studies  or  studies  involving 
systematic  changes  in  geometry  are  apt  to  have  limited 
success.   This  is  because  in-core  storage  limitations  do  not 
permit  subdividing  the  body  into  sufficiently  small  parts 
to  permit  obtaining  really  accurate  results.   Accordingly, 
one  of  the  recommendations  made  herein  is  that  attention  be 
focussed  on  modifying  this  program  so  as  to  permit  efficient 
out-of-core  storage  and  manipulation.   The  present  study  does 
demonstrate  that  the  procedural  sequence  adopted  and  imple- 
mented herein  leads  to  reliable  evaluations  if  the  geometry 
of  the  vibrating  object  is  sufficiently  simple  that  it  may 
be  represented  by  a  reasonably  small  number  of  finite 
el ements . 

Although  the  immediately  following  sections  (II,  III,  IV, 
V,  and  VI)  discuss  theory,  sources,  and  details  of  the  program 
developed  herein,  the  principal  thrust  and  value  of  this 


thesis  is  felt  to  be  the  detailed  instruction  provided  for 
the  employment  of  the  program  NAT! -3D  by  a  new  user  for  the 
actual  evaluation  of  the  fundamental  frequency  of  a  three- 
dimensional  elastic  body  in  which  he  may  be  interested. 


1 1 .   THE  FINITE  ELEMENT  METHOD 

In  recent  years  a  procedure,  called  the  Finite  Element 
Method,  has  been  developed  and  widely  used  for  a  variety  of 
engineering  problems  involving  the  analysis  of  bodies  having 
complicated  shapes.   In  the  analysis  of  the  stresses  and 
deformations  of  an  elastic  body,  using  this  method  (which  is 
generally  abbreviated  as  FEM),  the  body  is  subdivided  into 
a  (large)  number  of  geometrical  subdivisions  called  finite 
elements.   The  adjoining  surfaces  of  these  elements  are 
delineated  by  points  called  nodes.   In  a  three-dimensional 
problem  each  node  may  have  three  displacements.   Displace- 
ments, and,  by  a  process  of  differentiation,  strains  at 
points  interior  to  such  a  finite  element  are  expressable  in 
terms  of  nodal  displacements  by  use  of  so-called  shape 
functions  which,  roughly  speaking,  are  polynomials  in  the 
nodal  displacements.   The  laws  of  elasticity  are  used  to 
relate  nodal  displacements  of  a  finite  element  to  the  forces 
applied  at  its  nodes  and  the  relationship  is  conveniently 
expressed  in  the  form  of  an  element  stiffness  matrix. 

An  assembly  algorithm  is  used  to  combine  such  element 
stiffness  matrices  into  a  so-called  global  siffness  matrix, 
for  the  original  entire  body.   This  global  siffness  matrix, 
hereinafter  designated  as  [BGK],  relates  nodal  displacements 
for  the  entire  body  to  nodal  forces  which  may  be  applied  to 
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it.   If  there  are  a  large  number  of  finite  elements  in  the 
body  and  if  each  of  them  has  itself  a  substantial  number  of 
nodes,  the  total  number  of  displacements,  three  for  each  node, 
may  be  very    large. 

In  elastostatic  problems,  the  nodal  forces,  i.e.,  the 
force  components  applied  to  the  nodes,  are  arranged  in  the 
form  of  a  vector  and  the  unknown  nodal  displacement  components 
are  similarly  arranged,  so  that  the  relationship  can  be 
expressed  in  the  matrix  equation 

[BGK]  {v}  =  {f} 

If  the  size  of  the  vectors  {v}  and  {f}  were  not  great,  any 
of  a  number  of  procedures  could  be  used  to  solve  this  matrix 
equation  for  the  unknown  vector  {v}.   However,  in  practice, 
the  size  usually  is  quite  large,  and  severe  problems  of 
storage  in  an  electronic  computer  have  resulted  in  the 
development  of  quite  subtle  procedures  for  storing  the 
elements  of  [BGK]  and  of  {f}. 

The  most  significant  features  of  the  matrix  [BGK]  are 
that  it  is  positive  definite,  symmetric,  and  of  banded 
structure.   The  last  means  that  there  is  an  upper  right  and 
a  lower  left  triangle  which  consists  entirely  of  zero  elements 
which  need  not  actually  be  stored  if  their  presence  is  other- 
wise suitably  noted.   These  properties  permit  storing  the 
[BGK]  matrix  in  far  less  space  than  if  it  were  to  be  stored 
in  the  form  usually  used  to  exhibit  square  matrices. 
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Principle  features  of  any  FEM  computer  program  are  the 
generation  and  compact  storage  of  the  elements  of  [BGK],  and 
a  so-called  "equation  solver"  which  efficiently  solves  the 
linear  algebraic  problem  indicated  by  the  equation  above. 
Software  to  accomplish  this  has  been  developed  and  tested 
and  is  readily  available.   One  of  the  problems  facing  a  person 
who  undertakes  to  solve  an  engineering  problem  by  FEM  is  to 
select  among  several  competing  software  programs.   TRISOP 
[Ref.  3]  and  FIELD  [Ref.  4]  are  the  names  of  programs 
developed  by  Professor  G.  Cantin,  LCDR  E.A.  Leonidas  and 
LCDR  G.T.  Lew.   These  programs  have  been  successfully 
employed  at  NPS  for  a  variety  of  investigations.   It  was 
decided  to  adapt  and  employ  appropriate  sections  of  these 
programs  for  the  purposes  of  this  thesis. 

TRISOP  and  similar  programs  require  a  massive  input  to 
describe  the  geometry  of  the  body  and  the  manner  in  which  it 
is  subdivided,  to  introduce  the  physical  properties  of  the 
materials,  to  represent  the  manner  in  which  the  body  is 
constrained,  and  to  account  for  the  applied  loadings.   For 
example,  in  a  typical  problem  treated  in  the  course  of  this 
thesis  investigation  250  data  cards  would  be  required  as 
input.   Since  the  labor  of  generating  the  data  for  and 
actually  punching  these  cards  is  so  great  and  since  the 
likelihood  of  making  at  least  one  crippling  error  is  almost 
certain,  software  programs  called  "mesh  generators"  have 
been  developed  which  accept  a  minimum  of  input  data  and 


12 


instructions  and  which  produce  data  cards  suitable  for 
subsequent  use  in  an  FEM  program.   Such  a  mesh  generator, 
called  TRIMEG  [Ref.  2],  developed  by  LT  J.R.  Adamek,  for  use 
with  TRISOP  has  been  used  and  tested  at  NPS,  and  has  provided 
a  very  useful,  almost  essential  adjunct  to  the  FEM  program 
used  for  this  thesis. 

TRISOP  does  not  have  a  dynamic  capability.   As  will  be 
indicated  in  the  next  section  hereof,  one  of  the  important 
features  of  a  dynamics  FEM  program  for  use  in  determining 
natural  frequencies  is  an  appropriate  mass  matrix.   Without 
going  into  detail,  it  may  be  stated  that  the  theory  of  the 
FEM  provides  for  the  construction  of  so-called  consistent 
mass  matrices  for  the  individual  finite  elements  and  their 
assembly  into  a  global  mass  mass  matrix  [BGM].   The  computa- 
tional details  are  almost  exactly  the  same  as  for  the 
construction  of  element  stiffness  matrices  and  their  assembly 
into  [BGK].   The  global  mass  matrix  has  essentially  the  same 
properties  as  the  global  stiffness  matrix  and  requires  the 
same  storage. 

Much  thought  was  given  to  the  problems  of  storing  such 
a  consistent  mass  matrix.   However,  studies  by  others 
[Ref.  7]  have  indicated  that  when  evaluating  the  fundamental 
harmonic  frequency  of  free  vibrations  of  an  elastic  object, 
use  of  a  lumped-mass  form  of  the  mass  matrix  gives  results 
as  good  as  or  better  than  those  obtainable  by  using  a  con- 
sistent mass  matrix.   The  reasons  for  this  are  not  clear, 
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but  the  advantages  are.   A  lumped-mass  mass  matrix  is 
diagonal  in  form  and  can  be  stored  as  a  vector;  briefly, 
using  the  lumped-mass  form  results  in  a  dynamics  problem 
requiring  very  little  more  storage  than  an  elastostatic 
probl em  requi  res . 

Accordingly  it  was  decided  to  employ  the  lumped-mass 
form  and  the  question  arose  of  how  optimally  to  obtain  such 
a  mass  matrix.   The  procedure  which  finally  evolved  was  to 
construct  the  element  consistent  mass  matrix  and  to  diago- 
nalize  it  before  assembly  into  the  global  mass  matrix. 

The  word  diagnalize  in  this  context  is  not  intended  to 
mean  the  mathematical  process  of  di agonal i zati on  by  means  of 
solving  the  complete  eigenvalue  problem  for  the  matrix  M. 
Instead,  it  merely  means  replacing  the  consistent  (element) 
mass  matrix  by  a  diagonal  matrix  corresponding  to  an 
appropriate  division  of  the  mass  of  the  element  into  submasses 
concentrated  at  the  nodal  points. 

The  assembly  was  done  in  the  core  storage  space  reserved 
for  assembling  (and  later  operating  upon)  the  global  stiffness 
matrix,  and  the  assembly  algorithm  is  exactly  the  same  for 
both  matrices.   The  assembled  mass  matrix  is  diagonal  and 
its  nontrivial  elements  are  quickly  stored  elsewhere  in  the 
form  of  a  vector. 

The  details  of  di agonal i z i ng  the  element  mass  matrix 
deserve  a  brief  description,  explanation,  and  plausible 
justification.   Aside  from  the  fact  that  it  is  obviously 
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desirable  to  employ  a  procedure  which  involves  a  simple 
algorithm  and  results  in  the  same  diagonalized  element  mass 
matrix  regardless  of  renumbering  of  nodes,  a  reasonable  and 
plausible  distribution  of  the  total  element  mass  into  mass 
points  at  the  nodes  is  also  required.   This  matter  is 
discussed  in  somewhat  greater  detail  on  page  29. 

Accordingly  appropriate  segments  of  the  programs  TRISOP 
and  FIELD  were  selected  and  modified  so  as  to  incorporate 
the  construction  of  a  lumped-mass  mass  matrix  as  outlined 
in  the  preceding  paragraphs  and  so  as  also  to  incorporate 
the  operations  of  LDLT  decomposition,  matrix  (Stodola) 
iteration,  convergence  testing,  and  Rayleigh  quotient 
evaluation  discussed  in  the  following  section  hereof.   The 
resulting  program  is  called  NAT1-3D. 
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III.   METHOD  OF  SOLUTION 

With  the  use  of  the  [K]  and  [M]  matrices  provided  by 
the  finite  element  formation,  as  discussed  previously,  the 
free  vibration  problem  can  be  formulated  as 


[M]  {v}  +  [K]  {v}  =  {0} 


(1) 


The  remarks  which  follow  immediately  represent  standard 

procedure  in  vibration  theory  and  are  intended  merely  to 
recall  to  the  reader  the  operational  sequence  which  leads  to 

the  desired  result.   We  are  concerned  with  the  lowest 

frequency  harmonic  vibrations  which  the  system  may  execute. 

By  assuming  that  the  vector  {v  }  of  time-varying  displacements 
has  the  form 

{v}  =  {v*}cos(wt) 

(where  {v*}  is  a  vector  of  constant  quantities  representing 
the  amplitudes  of  the  motions  of  the  nodal  points)  and 
substituting  this  equation  into  equation  (1)  above,  the 
cosine  term  may  be  cancelled  out,  and,  upon  dispensing  with 
the  asterisk  and  henceforth  regarding  {v}  as  a  vector  of 
constants,  we  arrive  at  the  form 

[K]  {v}  =  032  [M]  {v} 
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It  may  be  shown  that  if  an  arbitrary  vector  {v.}  is 
introduced  on  the  right  side  and  this  equation  is  solved 
for  a  new  vector,  i.e.,  if  we  solve  for  {v.+,}  in  the 
equati  on 

[K]  Cvi+1}  =  u2[M]  {v.} 

then  {v. j,}  is  a  better  approximation  to  the  theoretically 
l  + 1 

correct  vector  for  the  lowest  frequency  than  was  {v.}. 
Furthermore,  if  the  vectors  are  suitably  and  similarly 
normalized,  the  scalar  constant  of  normalization  approaches 
more  and  more  closely  to  the  correct  value  of  the  square  of 
the  lowest  frequency.   This  process  of  matrix  iteration  is 
frequently  called  Stodol a-i terati on  after  the  Swiss  engineer 
who  introduced  the  method. 

Although  the  scalar  constant  of  normalization  does 

provide  an  increasingly  good  evaluation  of  the  desired  value 

2 

of  u-j  ,  at  any  stage  of  the  iteration  a  more  accurate  estimate 

may  be  obtained  from  the  Rayleigh  quotient 


2  ~ 


T 


T 


1 


=  Q  =  {v}'  [k]  {v}  /  {v}'  [M]  {v}. 


Performing  this  operation  after  a  reasonable  number  of 

iterations  provides  the  best  currently  available  value  for 

2 

W-,   and,  presuming  that  the  value  so  obtained  is  "close"  to 

the  values  previously  obtained  during  the  Stodola  iteration, 
provides  a  check  on  the  accuracy  and  the  validity  of  the 
entire  computational  process. 
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With  systems  having  only  a  few  degrees  of  freedom  it 
is  customary  to  solve  for  {v.  +  -.}  by  obtaining  the  inverse  of 
the  matrix  [K].   However,  the  desirable  property  of  the 
[K]  matrix  as  obtained  by  the  finite  element  method,  namely 
its  banded  character  which  permits  compact  storage,  is  not 
possessed  by  the  inverse  [K]~  ;  furthermore  the  operation 
of  obtaining  an  inverse  is  computationally  extravagant. 
Accordingly,  it  is  very    advantageous  to  represent  [K]  in  the 
so-called  LDLT  form,  i.e.,  to  decompose  [K]  into  the  matrix 
product 

[K]  =  [L]  [D]  [L]T 

where  [L]  is  a  lower  triangular  matrix  having  unit  diagonal 
elements  (which  need  not  be  stored)  and  [D]  is  a  diagonal 
matrix.   This  decomposition  is  a  standard  operation  which  is 
not  difficult  to  perform  and  the  matrices  [L]  and  [D]  may 
be  stored  in  the  space  occupied  by  [K]. 

The  normalization  which  is  used  is  the  so-called 
Euclidean  normalization  for  which  the  unnormalized  vector 
{v}  is  divided  throughout  by  the  square  root  of  the  sum  of 
the  square  of  its  elements. 

The  iteration  proceeds  as  follows.   Having  the  normalized 
vector  {v.}  (either  as  a  result  of  having  made  a  previous 
iteration  or,  the  first  time  through,  as  the  result  of  a 
guess  or  enlightened  estimate)  the  vector 


(g)  =  [M]  [v.] 
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is  obtained  simply  by  matrix  multiplication.   rhen  the 
vector  {h}  is  determined  from  the  equation 

[K]  fh}  -  (g) 

by  a  forward  and  backward  substitution  process  making  use 

of  the  desirable  properties  of  the  LDL7  decomposition.   Then 
{h}  is  normalized  so  as  to  produce  the  ne/t.  iterant  '7-,|}> 

the  reciprocal  of  the  constant  of  normalization  being  an 

2 
approximation  to  co-i  .   This  process  is  repeated  until 

convergence  is  satisfactory.   In  the  computer  Implementation 

of  this  procedure,  described  later,  convergence  was  assumed 

to  be  satisfactory  when  two  successive  such  evaluations  of 

7  -  fi 

oj ,  differed  by  less  than  10   .   At.  this  [joint,  the  Raylelgh 

quotient  was  evaluated  in  the  form 

Q  =  ([L]T  {v})T  [D]  LL]Tfv}  /  f  v  >T  [M]fvj 

Although  the  procedure  described  above  could  be 
employed  with  nondiagonal  [M],  the  fact  that  it  was  dec  id 
to  di  agonal  ize  [M]  in  order  to  minimize  corf:  storage 
requirements  also  made  the  computations  involving  [M]  very 
simple  and  economical. 
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I V .   PROGRAM  ARRANGEMENT 

To  facilitate  the  adaptation  of  existing  mesh  generating 
schemes  and  to  incorporate  the  strong  points  of  existing 
software  the  general  arrangement  of  the  program  is  modular 
in  nature.   This  also  allows  future  changes  to  be  accomplished 
more  easily  if  desired. 

The  actual  solution  of  a  problem  requires  two  basic 
steps:  mesh  generation  and  frequency  evaluation. 

A.    MESH  GENERATION 

Mesh  generation  is  the  process  in  which  the  object  to 
be  analyzed  is  divided  into  an  assemblage  of  small  elements. 
It  is  a  laborious  and  tedious  task  that  should  not  be  rushed 
into.   Intuition  guided  by  an  intelligent  study  of  the 
geometry  of  the  object  is  often  the  best  approach  to  use  in 
choosing  a  mesh. 

TRIMEG  (Ref.  1,  which  is  discussed  in  greater  detail 
in  Appendix  C)  is  a  computer  program  that  performs  mesh 
generation  for  three  dimensional  objects.   The  program  auto- 
matically discretizes  a  given  geometry  and  upon  command  will 
punch  data  decks  containing  nodal  coordinate  and  element 
connectivity  information.   The  solution  routine  presented  in 
this  thesis  was  written  to  accept  data  in  the  format  of 
TRIMEG's  output.   The  use  of  TRIMEG  is  not  a  necessity  and  any 
source  may  be  used  for  the  input  data.   However,  the  use  of 
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TR I  MEG  to  generate  data  cards  for  use  with  NAT1-3D  is 
strongly  recommended  to  save  effort  and  to  reduce  the  like- 
lihood of  error.   The  reader  is  advised  to  study  and 
understand  TRIMEG;  Ref .  1 . 

B.    FREQUENCY  EVALUATION 

1 .    Common  Statements 

Throughout  the  development  of  this  routine  a  prime 
objective  was  to  make  it  as  simple  as  possible  to  use. 
Simplicity  can  however  limit  flexibility.   To  enable  the 
user  to  avoid  requesting  excessive  core  space  by  adjusting 
the  core  requirements  to  meet  the  needs  of  his  particular 
problem,  the  large  data  blocks  are  placed  in  easily  changed 
COMMON  storage  areas.   This  adjustment  allows  the  use  of 
more  nodes  for  problems  with  smaller  band  widths. 

There  are  four  labeled  COMMON  statements  which  may 
be  used  to  modify  core  storage  requested  by  NAT1-3D;  this  is 
done  to  avoid  requesting  more  core  than  necessary  and  thus 
results  in  obtaining  optimum  job  priority.   The  locations 
of  the  cards  and  size  of  the  arrays  are  indicated  in  the 
listing  in  Appendix  B.   The  following  is  a  tabulation  of  the 
statement  names  and  the  subroutines  in  which  they  appear. 

BIGS  MAIN,  MERGE,  QUAD,  FORM,  FASTEN,  LTLD,  SOLV 

ENTER  INPUT,  MERGE 

ITER  MAIN,  NORM 

SOL  MAIN,  SOLV 
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In  addition  to  modifying  the  COMMON  statements  the 
dimension  of  the  vector  FMODE  which  appears  in  the  main 
routine  should  be    calculated  and  introduced  in  the  appropriate 
DIMENSION  statement. 
2.    Input 

The  required  data  consists  of  six  basic  card 
groupings.   Their  specific  formats  are  listed  in  Appendix  B. 

a.  TITLE  (one  card) 

Any  desired  problem  title  may  be  used. 

b.  PROBLEM  DATA  (one  card) 

The  card  contains  the  basic  problem  parameters. 
Here  the  user  chooses  the  type  of  shape  functions  that  are 
desired  and  indicates  the  number  of  Gauss  points  he  wishes 
to  use  in  their  integration. 

c.  MATERIAL  PROPERTIES  (one  card  per  material) 
These  cards  contain  the  elastic  property  data 

for  the  isotropic  materials  used  in  the  object.   The  number 
used  to  represent  a  material  is  arbitrary  but  it  must  be  the 
same  numbering  scheme  used  in  the  mesh  generation. 

d.  NODE  CARDS  (one  card  per  node) 

These  cards,  which  may  be  punched  by  TRIMEG, 
contain  the  node  numbers  and  their  global  coordinates. 

e.  CONNECTIVITY  MATRIX  (two  cards  per  element) 
These  cards,  which  may  be  punched  by  TRIMEG, 

contain  the  detailed  arrangement  of  all  the  contributions  of 
each  element  to  the  nodal  displacements.   As  the  name  implies. 
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this  matrix  connects  the  elements  and  serves  as  the  control- 
ling agent  which  distributes  the  contribution  of  each  element 
to  the  global  stiffness  and  mass  matrices. 

f.    BOUNDARY  CONDITION  CARDS  (one  card  per 
boundary  node) 

The  user  may  impose  boundary  conditions  by 

restraining  the  motion  of  any  node  in  any  coordinate 

direction.   The  boundary  condition  cards  indicate  the  node 

and  the  coordinate  direction  to  be  fixed. 
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V.   EXAMPLE  PROBLEMS 

The  method  usually  used  to  test  computer  programs  is 
to  solve  a  problem  that  has  a  known  exact  solution  to  serve 
as  a  comparison  gauge  for  the  computer  solution.   This 
presents  a  difficulty  as  the  author  is  unaware  of  a  truly 
three  dimensional  classic  solution  for  natural  frequency 
determination.   Problems  A  and  B,  below  may  seem  overly 
simple  but  they  do  have  known  approximate  solutions  and 
serve  to  establish  a  measure  of  confidence  in  the  integrity 
of  the  routine. 

The  approximate  solutions  are  based  on  Euler  beam  theory 
or  on  the  more  elaborate  Timoshenko  beam  theory,  both  of 
which  involve  replacing  the  equations  of  three  dimensional 
elasticity  (theoretically  required  to  describe  the  problem) 
by  greatly  simplified  equations. 

A.    CANTILEVER  BEAM 

The  cantilever  beam  is  a  basic  design  structure  which 
may  vary  in  form  from  the  smallest  man-made  gear  tooth  to 
the  largest  trees  produced  by  mother  nature.   The  difficulty 
of  analysis  also  may  vary  as  drastically.   The  inclusion  of 
effects  such  as  shear  and  rotary  inertia  significantly 
complicates  the  frequency  equation. 


24 


The  specific  structure  used  in  this  problem  was  a  beam 

six  inches  long  having  a  uniform  one  inch  square  cross 

section.   Subdivision  is  indicated  in  Figure  1.  The  results 
were  as  follows: 


Euler  Beam  Theory 
NAT1-3D  Results 
Ratio  NAT1-3D/Euler 


a)  =  5700  rad/sec 
w  -    5664  rad/sec 


0.994 


This  agrees  with  Timoshenko  theory  which  predicts  a 
frequency  correction  for  shear  and  rotary  inertia  of  nearly 
unity  for  this  geometry  [Ref.  8].   Figure  1  shows  the 
convergence  of  the  solution  as  the  number  of  nodes  is 
increased,  i.e.,  the  effect  of  finer  subdivision.   Note 
however,  that  the  finer  discretization  involves  only  the 
"Z"  direction.   It  would  be  of  interest  also  to  investigate 
finer  subdivisions  in  the  "X"  and/or  "Y"  direction.   However, 
this  would  involve  substantial  increases  in  bandwidth. 

B.    SIMPLY  SUPPORTED  BEAM 

The  effects  of  shear  deflection  and  rotary  inertia  are 
more  pronounced  in  a  simply  supported  beam  having  the  same 
geometry  as  the  previous  cantilever. 

Euler  Beam  Theory        a  =  16,000  rad/sec 
NAT1-3D  to  =  14,250  rad/sec 

Timoshenko  Beam  Theory   u>  =  14,000  rad/sec 
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The  Timoshenko  solution  is  based  on  an  assumed  shear 
constant  value  of  0.6.   This  value  is  the  approximate 
middle  of  the  accepted  range  of  values  [Ref.  8]. 

C.    CANTILEVER  BEAM  WITH  COMPLIANT  FOUNDATION 

Figure  2  represents  a  mesh  which  could  be  used  to 
investigate  the  effects  of  root  fixity  on  the  natural 
frequency  of  a  cantilever  beam  having  imperfect  root  fixity 
The  boundary  conditions  would  be  applied  to  the  outer  wall 
surfaces  thus  allowing  the  root  to  deflect  naturally.   The 
material  properties  of  the  wall  and  the  beam  can  also  be 
varied  to  measure  the  effects  of  dissimilar  materials.   A 
final  solution  was  not  obtainable  due  to  the  large  core 
requirements  arising  from  even  the  rather  coarse  discretiza 
tion  indicated  in  Figure  2. 
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VI.   DISCUSSION 

A.    MASS  MATRIX 

As  noted  in  chapter  II  a  lumped  diagonal  mass  matrix 
was  chosen  for  use  in  the  solution  routine.   The  choice  was 
based  upon  the  reduced  core  requirements  and  computational 
advantages  inherent  to  the  diagonal  matrix.   A  decision  was 
made  to  form  the  lumped  matrix  from  the  consistent  matrix 
because  the  numerical  evaluation  of  the  consistent  mass 
matrix  parallels  the  evaluation  of  the  stiffness  matrix 
allowing  the  two  formulations  to  be  coded  simultaneously. 
Lumping  is  accomplished  by  replacing  the  diagonal  entry  in 
each  row  of  the  consistent  matrix  with  the  sum  of  the  entries 
in  that  row  and  zeroing  the  non-diagonal  elements.   The 
result  could  be  termed  a  consistent,  lumped  matrix. 

Reference  7,  previously  cited,  indicate  that  this  form 
of  mass  matrix  yields  excellent  results  for  the  lowest 
frequency.   It  is  of  some  interest  at  this  point  to  examine 
the  point-mass  distribution  which  results  from  this  process 
in  the  case  of  an  element  in  the  form  of  a  rectangular 
parallelepiped,  a  case  for  which  physical  intuition  provides 
some  indication  of  whether  or  not  the  substitution  is 
reasonabl e . 
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Figure  3  shows  the  resultant  mass  distribution  for  both 
linear  and  quadratic  elements.   The  distribution  for  the 
linear  element  appears  logical.   The  quadratic  element 
however  is  not  at  all  obvious.   The  allocation  of  negative 
mass  at  the  element  corners  would  seem  to  conflict  with 
common  sense.   These  results  are  however,  not  without 
precedent.   Negative  corner  values  result  when  imposing 
uniform  body  forces  on  rectangular  elements  in  a  pattern  that 
is  nearly  identical  to  that  shown  in  Figure  1  of  Ref.  6 
(p.  168). 

B.    SHAPE  FUNCTIONS  USED 

In  its  earlier  versions  NAT1-3D  accommodated  quadratic 
shape  functions  only.   This  was  a  compromise  between  the 
limitations  of  linear  functions  in  representing  complex 
geometry  and  the  large  core  size  requirements  of  cubic 
elements.   A  later  inclusion  of  linear  shape  functions 
resulted  in  the  observation  that  they  were  "too  stiff."   The 
answers  they  gave  for  the  fundamental  frequency  were  high 
when  the  same  number  of  nodes  and  same  nodal  arrangement  was 
used  that  corresponded  to  an  accurate  solution  based  on 
quadratic  functions. 

This  is  to  be  expected  since,  for  the  same  subdivision 
into  elements,  the  shape  functions  for  quadratic  elements 
permit  a  closer  representation  of  the  actual  deformation 
pattern  than  is  permitted  by  using  linear  elements.   In 
effect  this  corresponds  to  both  types  of  elements  implicitly 
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Mass  ,  Rectangul ar 
Parallelepiped, 
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Distribution  of  Mass, 
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Quadratic  Functions 
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Figure  3.   Allocation  of  Lumped-Mass  Mass  Matrix 
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introducing  fictitious  and  undesirable  restraints,  which 
increase  the  calculated  frequency  above  its  correct  value; 
however,  these  restraints  are  greater  and  have  a  more  serious 
effect  in  the  case  of  linear  elements  than  they  do  in  the 
case  of  quadratic  elements. 


C.    NUMBER  OF  GAUSS  POINTS 

If  N  represents  the  number  of  Gauss  points,  the  Gaussian 
technique  used  in  this  routine  will  exactly  integrate  a 
polynomial  of  degree  (2N-1).   It  has  been  noted  however,  in 
calculations  employing  NAT1-3D  that  the  minimum  number  of 
Gauss  points  which  can  successfully  be  used  with  quadratic 
elements  is  three.   When  fewer  than  three  Gauss  points  are 
used  it  is  found  that  accuracy  of  the  frequency  evaluation 
is  seriously  decreased.   It  is  conjectured  that  this  behavior 
is  related  to  the  method  of  introducing  the  boundary  condi- 
tions (by  multiplying  selected  diagonal  elements  of  [BGK] 
by  a  large  factor,  "WELD").   The  calculations  which  showed 
the  undesirable  behavior  had  WELD  =10   .   It  would  be 
interesting  to  look  into  this  matter  more  fully  by  reducing 
the  value  of  WELD,  by  considering  alternate  methods  of 
introducing  boundary  fixity,  and  in  other  ways.   Briefly,  we 
do  not  presently  have  a  satisfactory  understanding  of  the 
effect  of  modifying  the  number  of  Gauss  points  except  to 
observe  that  using  NGP  =  2  gave  unsatisfactory  results,  while 
using  NGP  >  3  gave  results  which  were  felt  to  be  no  more 
accurate  or  reliable  than  those  obtained  by  using  NGP  =  3. 
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D.  CONVERGENCE 

Our  experience  has  shown  that  the  Stodola  iteration 
converges  with  amazing  rapidity,  three  to  four  iterations 
usually  being  sufficient.   This  can  be  attributed  to  the 
fact  that  the  structures  considered  had  lower  frequencies 
which  were  widely  separated  from  the  second  frequency.   It 
is  expected  that  objects  whose  lower  frequencies  are  not 
separated  will  require  more  iterations  to  converge. 

It  should  be  noted  that  the  "initial  guess,"  being  a 
unit  vector,  contains  elements  of  both  transverse  and  axial 
vibrations.   If  the  lowest  mode  were  an  axial  mode  it  would 
be  reflected  in  the  corresponding  elements  of  the  stiffness 
matrix  and  the  iteration  should  converge  to  the  actual 
lowest  frequency  be  it  an  axial  or  transverse  mode.   Correct 
solutions  have  been  obtained  for  transverse  vibration  using 
a  pure  axial  mode  as  an  initial  guess. 

E.  PROGRAM  LIMITATIONS 

It  could  be  considered  overkill  to  use  this  routine  on 
a  truly  two  dimensional  plane  strain  or  plane  stress  problem 
During  the  development  of  this  routine  a  2-D  version  was 
written  but  it  required  exactly  the  same  work  on  the  user's 
part  in  the  form  of  data  preparation  and  further  required 
the  user  to  recognize  that  his  object  could  be  modeled  as  a 
plane  stress/strain  case.   Upon  considering  that  an  original 
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goal  was  to  limit  the  assumptions  and  idealizations  that  the 
user  has  to  make,  only  the  three  dimensional  routine  is 
being  presented. 

In  its  present  format  the  program  is  limited  by  core 
size.   This  problem  arises  when  dealing  with  complex  geome- 
tries which  require  several  layers  of  elements  in  more  than 
one  direction.   Multidirectional  layering  results  in  a 
dramatically  rapid  growth  in  the  band  width  of  the  stiffness 
matrix  which  is  the  major  core  consuming  element  of  the 
entire  routine. 
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VII.   RECOMMENDATIONS 

A.  ELIMINATE  STORAGE  LIMITATIONS 

The  first  and  most  obvious  recommendation  for  further 
improvement  of  the  routine  is  to  remove  the  core  limitations 
imposed  on  the  working  area.   This  can  be  accomplished  by 
the  use  of  direct  access  devices  and  a  block  solving  routine 
which  secti onal i zes  the  matrices  and  works  in  core  with  only 
one  section  at  a  time. 

It  should  be  noted  that  the  implementation  of  these 
modifications  will  be  a  difficult  task.   Substantial  coding 
changes  will  be  required  and  the  difficulties  encountered 
are  anticipated  to  be  orders  of  magnitude  greater  than  those 
associated  with  the  in-core  routine. 

B.  PARAMETER  STUDY 

It  is  normally  a  part  of  verifying  and  testing  a  FEM 
program  to  perform  parametric  studies  on  several  problems 
having  known  exact  or  approximate  results.   Time  did  not 
permit  doing  such  an  extensive  investigation  with  NAT1-3D 
and  it  is  recommended  that  appropriate  attention  be  given  to 
a  systematic  study  of  the  effects  of  varying  the  values  of 
NGP,  WELD,  and  the  degree  of  the  shape  functions  and  of 
varying  the  fineness  of  the  ultimate  discretization  in  each 
of  the  three  coordinate  directions.   The  literature  on  FEM 
describes  many  such  studies  and  can  be  used  to  guide  such  an 
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investigation.   However,  it  may  be  stated  that  such  previous 
investigations  have  not  led  to  any  significant  understanding 
of  the  effects  of  varying  parameters  in  general.   Apparently, 
as  of  yet,  it  is  important  to  carry  out  such  an  investigation 
for  each  new  FEM  code. 
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APPENDIX  A 
PROGRAM  DESCRIPTION 

A  complete  disection  and  detailed  description  of  the 
program  is  not  practical.   This  appendix  is  intended  to 
contain  information  which  will  assist  both  a  reader  who 
desires  to  pursue  a  more  in-depth  study  of  the  problem  and 
a  reader  who  wishes  only  to  make  an  intelligent  casual  use 
of  the  program. 

A.    TERMINOLOGY 

The  following  is  a  listing  and  brief  description  of 
important  constant,  vector ,  and  matrix  names.  Functions, 
definitions  and  dimensions  are  given  if  applicable. 

1 .    Constants 

FACT  Value  of  the  Euclidean  norm 

of  the  mode  vector,  used  to 
determine  iteration  convergence 

NBAND  Half  band  width  of  the  stiff- 

ness matrix 

NEL  Number  of  elements 

NEQ  Number  of  equations, 

NEQ  =  3  x  NNODE 

NGP  Number  of  Gauss  points  used 

in  the  Gaussian  quadrature  . 
formul a 

NMAT  Number  of  materials 

NNBC  Number  of  bounded  nodes 

NNODE  Total  number  of  nodes 

NPEL  Number  of  nodes  per  element 
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2.    Vectors  (size) 

B 

EMODE  }     (NEQ) 
FMODE 

BGM 


Working  vectors  used  to  store 
the  mode  shape  during  iteration 
and  while  evaluating  the 
Rayleigh  quotient 


Di agonal i zed  mass  matrix  stored 
as  a  vector 


3.    Matrices  (dimensions) 


BGK(NEQ,NBAND) 
C0ARD(NN0DE,3) 

ENW(3*NPEL,3*NPEL) 
NC0HN(NEL,NPEL+1  ) 
NBC(30,4) 


Global  or  system  stiffness  matix 

Matrix  containing  global 
coordinates  of  nodes 

Element  stiffness  or  mass  matrix 

Connectivity  matrix 

Matrix  of  bounded  nodes  (i.e., 
nodes  on  which  boundary  condi- 
tions are  to  be  imposed)  and 
their  boundary  conditions 


B.    PROGRAM  ARRANGEMENT 

As  indicated  in  Figure  4  the  routine  is  of  modular 
construction.   This  format  allows  for  maximum  ultilization 
of  existing  software  and  facilitates  future  revisions. 

1.    MAIN 

The  first  segments  of  the  MAIN  program  initiate  the 
calculation  of  the  large  system  matrices.   Data  is  passed  to 
the  various  subroutines  via  common  statements.   The  iterative 
portion  of  MAIN  is  based  upon  the  decomposition  of  BGK  and 
the  repeated  solution  of  the  matrix  equation 
[BGK]  {v}  =  u  [BGM]  {v}.   Once  convergence  has  been  attained 
MAIN  controls  the  calculation  of  the  Rayleigh  quotient  which 
serves  as  both  a  check  and  a  refinement  of  the  frequency 
obtained  by  the  iteration. 
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ARRANGEMENT  OF  NAT1-3D 
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NORM 


RAYLEIGH  QUOTIENT 


■>  (STOP  on  empty  data  file) 


NAMES  in  boxes  are 
Subroutine  names 


Figure  4 
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2.  Subroutine  INPUT 

INPUT  reads  and  echo  checks  the  problem  data.   The 
most  common  source  of  errors  in  a  large  routine  is  a  simple 
typographical  error  in  the  data  deck.   INPUT  also  computes 
the  half  band  width  of  the  system  matrices. 

3.  Subroutine  MERGE  [Ref.  4] 

MERGE  combines  the  contributions  of  each  element 
and  distributes  them  to  the  corresponding  global  matrix, 
[BGK]  or  [BGM].   The  distribution  allocation  is  controlled 
by  the  connectivity  matrix  (NCONN).   MERGE  serves  in  a 
sequence  controlling  capacity  throughout  the  calculation  of 
the  system  stiffness  and  mass  matrices.   BGM  is  computed 
first  and  during  its  assembly  is  stored  in  the  region  set 
aside  for  BGK.   Calculation  of  the  weight  of  the  object 
provides  a  quick  check  on  the  mass  matrix  by  giving  the  user 
an  evaluation  which  he  may  compare  to  the  actual  weight  of 
the  object. 

4.  Subroutine  QUAD  [Ref.  4] 

QUAD  is  divided  into  two  basic  sections  which 
differ  in  the  form  of  the  shape  functions  used.   Computing 
the  mass  matrix  requires  evaluation  of  the  shape  functions 
while  computing  the  stiffness  matrix  requires  evaluation  of 
the  derivatives  of  the  shape  functions.   QUAD  selects  the 
weighting  factors  and  integration  points  that  correspond  to 
the  number  of  Gauss  points  desired  by  the  user.   It  then  calls 
the  appropriate  entry  of  Subroutine  FORM  which  evaluates  the 


40 


integrand  of  the  Gaussian  quadrature  formula.   Next  QUAD 
multiplies  the  result  by  the  Gaussian  weighting  factor  and 
computes  the  summation.   Di agonal i zati on  of  the  mass  matrix 
is  also  accomplished  in  QUAD. 

5.  Subroutine  ELASM  [Ref.  3] 

ELASM  forms  the  constitutive  elasticity  matrix  for 
each  element  based  on  its  material  which  is  assumed  to  be 
linear,  isotropic  and  elastic. 

6.  Subroutine  FORM  [Ref.  4] 

FORM  evaluates  the  product  of  the  shape  functions 
or  the  product  of  their  partial  derivatives  at  each  selected 
Gauss  point  in  the  local  coordinate  system.   Dual  entry 
points  are  used  to  distinguish  between  the  computations 
necessary  to  form  the  elemental  mass  matrix  and  the  elemental 
stiffness  matrix. 

7.  Subroutine  SHAPE  [Ref.  4] 

SHAPE  performs  the  actual  evaluation  of  the  shape 
functions  and  their  derivatives  with  respect  to  the  local 
coordinate  system.   Provisions  are  included  for  both  linear 
and  quadratic  shape  functions. 

8.  Subroutine  JACOB  [Ref.  4] 

JACOB  evaluates  the  Jacob i an  matrix  as  well  as  its 
inverse  and  its  determinant. 
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9.    Subroutine  FASTEN 

FASTEN  spatially  fixes  the  bounded  nodes  in  the 
desired  coordinate  direction  by  multiplying  the  appropriate 
diagonal  elements  of  [BGK]  by  10    so  as  effectively  to 
suppress  the  associated  deflection  component.   This  is  a 
commonly  used  method,  in  FEM  analysis,  of  introducing  fixity 
at  boundary  nodes. 

10.  Subroutine  LDLT  [Ref.  5] 

LDLT  decomposes  the  system  stiffness  matrix  [BGK] 
into  the  products  of  three  matrices:   a  lower  triangular 
matrix  [L]  having  unit  values  on  the  diagonal,  a  diagonal 
matrix  [D],  and  [L]   the  transpose  of  [L];  viz. 
[BGK]  =  [L][D][L]T.  In  the  process  [BGK]  is  destroyed  and 
replaced  by  the  elements  of  [D]  and  the  nontrivial  elements 
of  [L].   The  diagonal  matrix  is  stored  in  the  first  column 
and  the  unit  triangular  matrix  is  stored  in  the  remaining 
space . 

11.  Subroutine  SOLV  [Ref.  5] 

SOLV  performs  a  forward  substitution  followed  by 
a  back  substitution  to  solve  the  system  [BGK]  {X}  =  {B}. 
The  answers,  {X},  are  returned  in  vector  {B}. 

12.  Subroutine  NORM 

NORM  calculates  the  Euclidean  norm  of  the  mode 
vector.   This  norm  is  equal  to  the  square  root  of  the  sum 
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of  the  squares  of  the  element  values  in  the  vector.   Each 
element  in  the  vector  is  then  divided  by  the  norm  which 
results  in  normalizing  the  vector. 

C.    CORE  ESTIMATION 

The  mesh  generating  scheme  chosen  by  the  user  should 
inform  him  of  the  number  of  nodes  and  the  bandwidth  that 
have  resulted  from  the  discretization.   TRIMEG  clearly  lists 
this  data.   With  the  aid  of  the  following  equation  these 
figures  may  be  used  to  estimate  the  core  requirements  which 
result  from  the  mesh. 

"K"  required  .  110  +  24  x  NNODES^NBAND  +  4) 

The  resultant  answer  includes  the  necessary  buffers  and 
can  be  used  directly  as  the  region  specified  on  the  FORTCLG 
control  card. 
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APPENDIX  C 
MESH  GENERATION  USING  TRIMEG 

The  purpose  of  this  appendix  is  to  give  the  user  who 
wishes  to  use  TRIMEG  a  few  brief  but  very    important 
suggestions  that  could  save  time  and  reduce  frustration. 

A.  BACKGROUND 

TRIMEG  [Ref.  1]  in  an  automatic  mesh  generator  for 
use  with  three  dimensional  isoparametric  elements.   It 
is  both  quick  and  efficient  and  saves  the  user  a  lot  of 
tedious,  time  consuming  work  by  punching  the  bulk  of  the 
data  cards  required  to  solve  the  vibrational  problem. 

B.  COMMENTS 

When  using  TRIMEG  there  are  a  few  subtle  guides  which 
must  not  be  overlooked  in  order  to  obtain  the  most  efficient 
discreti  zati  on . 

1 .    Orientation 

Unless  the  user  has  a  reasonably  good  idea  of  the 
purposes  and  functions  of  TRIMEG  it  is  likely  that  he  will 
not  use  NAT1-3D  to  the  greatest  advantage.   He  is  likely  to 
orient  his  super  elements  or  provide  for  subsequent  sub- 
division of  super  elements  into  a  number  of  individual 
elements  in  a  disadvantageous  way  so  as  to  use  up  available 
storage  while  accommodating  a  smaller  total  number  of  elements 
than  would  have  been  possible  by  a  more  effective  arrangement. 
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Accordingly  he  is  urged  to  study  Ref.  2  before  undertaking 
a  problem.   In  particular,  referring  to  Figure  13,  p.  29, 
Ref.  2,  the  following  general  remarks  may  be  of  assistance. 

The  super  elements  all  have  the  same  orientation  given 
by  xi ,  eta,  and  zeta  directions.   A  super  element  is  divided 
into  ROW  parts  in  the  eta  direction,  into  COL  parts  in  the 
xi ,  direction,  and  into  SLICE  parts  in  the  zeta  direction. 
The  body  represented  by  the  super  elements  has,  let  us  say, 
a  maximum  of  nxi  super  elements  in  the  xi  direction,  a 
maximum  of  neta  super  elements  in  the  eta  direction,  and  a 
maximum  of  nzeta  super  elements  in  the  zeta  direction.   The 
orientation  of  the  super  elements  should  be  chosen  so  that 
the  numerical  products  (R0W)(nxi),  (C0L)(neta),  and  (SLICE) 
(nzeta)  are  in  ascending  order.   This  gives  minimum  bandwidth 
corresponding  to  the  total. number  of  elements  into  which  the 
body  is  finally  subdivided. 

2.    Boundary  Nodes 

For  configurations  having  curved  boundaries  it  is 
possible  and  advisable  to  use  32  boundary  nodes  to  describe 
the  super  element  even  when  linear  or  quadratic  shape 
functions  are  to  be  used.   The  use  of  the  additional  nodes 
allows  TRIMEG  to  more  accurately  follow  the  curvature  of  the 
boundary  of  the  object  as  it  subdivides  the  super  element. 
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LISTING    OF    DATA    DECK    FOR    TRIMEG 
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